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ABSTRACT 



This thesis investigates the behavior of Graphite- Epoxy 
composites subjected to fires as may occur on the decks of 
naval aircraft carriers. The analytical model consists of two 
parts; one for the determination of the temperature field 
within the composite due to a fire and the other for 
determining the stresses within the composite due to the 
temperature field. Both problems are provided one -dimensional 
finite element models. Appropriate failure criteria are 
incorporated to predict the survivability of composites in 
various fire environments. Parametric studies were performed 
with the results presented in both graphical and tabular form. 
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I. 



INTRODUCTION 



A. BACKGROUND 

Composite materials, due to their high strength to weight 
ratio, are becoming increasingly common as construction 
materials in today's modern combat aircraft. Through use of 
these materials the weight to strength ratio of the airframe 
has been greatly reduced allowing for an increase in the 
respective payload per unit weight of the structure. 

Several different Navy fixed wing aircraft such as the F- 
18, F-16, F-15 and F-14 all use composite materials in their 
construction. Since these aircraft are used on board naval 
aircraft carriers, it is quite conceivable that these 
structures could be exposed to flight deck, fires of varying 
intensities during their service life. It is important for 
the Navy to assess the survivability of composite aircraft 
exposed to fires. To this end, experimental and analytical 
investigations have been undertaken. 

In 1979 testing was performed at the Naval Weapons Center 
(China Lake) to determine the relative survivability of 
Graphite -Epoxy Composite materials and 7075 -T6 Aluminum Alloy 
when exposed to flight deck pool fires [Ref. 1] . These tests 
indicated the susceptibility to burn through for Graphite- 
Epoxy samples was approximately the same as for the comparable 
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7075 -T6 Aluminum Alloy samples tested. However unlike the 
aluminum, in some cases the composite material continued to 
smolder and/or burn after application of extinguishing agents. 
In an attempt to model this smoldering/burning of Graphite- 
Epoxy Composites, Vatikiotis in 1980 [Ref. 2] developed a 
computer model which he later refined in 1982 [Ref. 3] . This 
finite element program evaluates fiber diameters, fiber 
temperatures, air temperatures and oxygen levels throughout 
Graphite- Epoxy materials exposed to fire. The evaluations are 
with respect to time and location. 

Graphite -Epoxy material consists of high strength graphite 
fibers surrounded by an epoxy matrix. When exposed to a fire 
the epoxy matrix material burns off rather quickly at 400- 
500°F leaving only the graphite fibers supported by residual 
combustion products. These remaining high strength fibers 
develop thermally induced stresses due to the varying 
temperature profile across the material. If high enough 
stresses are developed, fibers may fail rendering the 
remaining material of no structural use. However, if a 
considerable amount of fibers do not fail the surviving 
material, although considerably weakened in the compression 
mode, could continue to be used in the tensile mode. 

B, OBJECTIVE OP THIS THESIS 

This thesis investigated thermally induced stress failure 
of Graphite- Epoxy material exposed to varying intensity fires. 
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The analysis was carried out in two parts. In part one, the 
computer model developed by Vatikiotis in 1982 was used to 
model fires of different intensities. This program evaluated 
the temperature profiles across a porous medium consisting of 
graphite fibers as time progressed. These results were 
further utilized in the second part of the analysis. Part two 
of the analysis consisted of developing a stress model and 
finite element computer program. This stress code evaluated 
the thermally induced stresses of fibers based on the 
temperature profiles determined from the combustion code. The 
fiber stresses were calculated at each time step; failure 
criteria applied and failed fibers removed from the stress 
analysis. Through this procedure, a quantitative evaluation 
of Graphite -Epoxy material failure with respect to time was 
obtained. 

C. THESIS OUTLINE 

Chapter II presents the physical and mathematical models. 

Chapter III develops the failure criteria used in 
determination of individual fiber failure. 

Chapter IV contains a description of the stress analysis 
program and it's interaction with the combustion program 
developed by Vatikiotis. 

Chapter V presents case studies with results presented in 
graphical form. 

Chapter VI contains the conclusions. 



3 



II. PROBLEM FORMULATION 



A. PHYSICAL MODEL 

The physical model assumed for this study is that of a 
Graphite -Epoxy composite plate exposed to fires of varying 
intensity as shown in Figure 1. This plate of composite 
material could have several different nominal thicknesses 
depending on the number of laminate layers or ply assumed. 
For this study, nominal thicknesses of 0.25, 0.5, 0.75 and 1.0 
inch are assumed, corresponding to approximately 35, 70, 105 
and 140 ply composites respectively. 
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As the composite plate heats up it will lose its epoxy- 
matrix material between 400 and 500°F. The matrix loss 
process is too complex to model exactly. However, for epoxy 
loss to occur, the epoxy matrix must undergo phase changes 
from solid to combustion products. These phase changes occur 
at a constant temperature throughout the material . 

This study assumes that no failure of graphite fibers 
occurs prior to the burn- off of the epoxy matrix material. The 
point at which the epoxy matrix has been reduced to residue 
will be chosen as the starting time for all case studies. 
Because of the phase changes which occur during the matrix 
burn off process described above, the initial temperature 
profile across the remaining Graphite fiber medium is assumed 
to be a constant 400°F at the starting time. 

Air flows across the upper surface of the composite plate. 
This surface air flow produces a differential pressure across 
the thickness of the plate. At the starting time for these 
studies, since the matrix material has been removed, this 
differential pressure will induce air flow perpendicular 
through the remaining composite plate. 

B. VATIKIOTIS COMBUSTION MODEL 

In order to obtain the temperature field required for 
stress analysis the one -dimensional finite element program 
developed by Vatikiotis was used. As seen in Figure 2 his 
combustion model assumes the composite material has been 
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exposed to heat sufficient to remove the epoxy matrix leaving 
only the graphite fibers surrounded by air filled voids. The 
heat supplied by the fire is implemented through application 
of a constant heat flux boundary condition at the lower, X=0, 
surface. Differential pressure across the plate thickness 
forces air flow through the material. Prior to fiber 
combustion the incoming air is considered to be of standard 
composition containing no fire combustion products. After 
combustion begins carbon dioxide and carbon monoxide also flow 
through the material. 
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Figure 2: The Vatikiotis combustion model. 



In his presentation Vatikiotis divides the material across 
its thickness into 30 equal thickness elements having a total 
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of 31 nodal points. He further develops a system of four 
transient one -dimensional equations consisting of an energy 



balance on the 


graphite 


material , 


an energy 


balance 


on 


the 


internal air. 


conservation of 


species on 


oxygen 


and a 


combination of 


Darcy' s 


law and 


continuity 


equation 


for 



internal air flow velocities. 

The two energy balances considered a differential volume 
of material with the following convention: 

Heat into ^ Heat _ Heat out of ^ Increased 

dv Generation dV Internal Energy 

The energy balance for the graphite fibers accounted for 
heat transfer into and out of the differential volume by three 
modes; conduction, convection and radiation, using first term 
Taylor series approximations. The heat generation term, 
resulting from burning of the graphite fibers, was modeled 
using an expression of Arrhenius type. 

The energy balance for the internal air considers only 
conduction and convection at the differential volume surfaces 
along with energy transport due to the flow of air through the 
volume. Also this equation does not contain a heat generation 
term. 

The oxygen mass balance started from the basic form of: 

Oxygen into ^ Oxygen out of ^ Oxygen ^ Oxygen 

dv dv Consumption Accumulation 

The oxygen into and out of the differential volume 
considers both molecular diffusion and convective transport 
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across the differential volume surfaces, again using a first 
term Taylor series approximations. The oxygen consumption was 
based on the graphite fibers (carbon) burning with the oxygen 
present, in a appropriately specified stoichiometric ratio. 
This reaction is assumed to produce carbon dioxide and carbon 
monoxide only. A relationship for the consumption of oxygen 
was developed through combination of the stoichiometric ratio 
with the Arrhenius reaction rate equation. 

The final equation of this system is the equation for 
internal air flow velocities. This equation was developed 
using a combination of Darcy's law for flow in a porus medium 
and the continuity equation. 

The above four equations form a set of coupled, nonlinear, 
transient, partial differential equations. These equations 
are solved for graphite temperature, air temperature, oxygen 
concentration, and internal air flow velocities with respect 
to time and position. Incorporated in this mathematical model 
are the variations of temperature dependent properties namely, 
graphite thermal conductivity, fluid (air) viscosity, and air 
density, along with parameters derived from them. These 
properties and parameters are recalculated at each time step 
of the transient solution. The Darcy's law-continuity 
equation was solved using a shooting method and the remaining 
three equations were solved using a Galerkin formulation of 
the finite element method. For completeness the above 
equations in their final form along with their associated 
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boundary conditions are contained in Appendix A. For more in 
depth review the reader is referred to Vatikiotis' original 
works . 

The main observations of Vatikiotis' earlier work are very 
obvious when his computer code is used. The differential 
pressure across the porous medium causes a flow of air to flow 
through medium. This induced air flow provides for both 
convective heat transfer of the fibers and a ready supply of 
oxygen for the combustion process. These two effects of the 
air flow are conflicting effects. First, the supply of oxygen 
increases the rate of combustion which increases the fiber 
temperature. Second, the convective heat transfer away from 
the hot fibers reduces the fiber temperature. A system 
comprised of a plate and a fire is defined by (a) the fire 
intensity and its time duration and (b) the plate subjected to 
the fire. A plate with known initial temperatures and oxygen 
concentrations, thickness, porosity, and pressure differential 
is subjected to a fire of specified intensity for a given time 
period. When the fire is removed the plate may either 
continue to burn (plate combustion) or the plate may return to 
ambient conditions (plate extinction) . Which process occurs 
depends on the state of the plate when the fire has been 
removed. If the temperature and air velocity fields exceed a 
critical state then the convective heat transfer from the 
fibers to the air is overtaken by the heat generated by the 
fiber combustion and the material will continue to burn. If 
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the critical state was not achieved, the convective heat 
transfer will predominate resulting in cooling and extinction. 

C. STRESS ANALYSIS MODEL 

The combustion model described previously determines the 
temperature distribution across the thickness of the medium. 
The stress analysis model uses the temperature profiles from 
the combustion code to determine the thermally induced 
stresses in the fiber composite plate. 

As seen in Figure 3 the stress model was developed 
assuming the graphite fibers were cylindrical rods running 
parallel to the plate surfaces. These fiber rods are modeled 
as rigidly fixed at one end and on the opposite end are 
modeled to move in unison with the other fibers. The 
expansion or contraction of fibers due to varying temperatures 
will cause forces to develop in the other fibers (See Figure 
4) . These resulting forces could cause excessive stresses to 
develop possibly causing failure. The final solution to this 
statically indeterminate problem required use of finite 
element techniques. 
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As seen in Figure 5 each element (fibers per ply) was 
modeled as a linear elastic spring with forces applied at each 
end (nodal point). The finite element "direct approach" [Ref. 
7] was applied to each individual element to obtain the 
element stiffness matrix and force vector as follows. All 
forces, displacements, stresses, strains and element 
characteristics are in the axial direction only. 




First referring to Figure 5, using summation of forces: 

5^F=F^ + F2 = 0 (2.1) 

F^ = -F 2 (2.1a) 

where : 

Fj^ = force applied at point 1 
F 2 = force applied at point 2 
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Assuming a linear spring model yields: 



F2=k{u^-u^) ( 2 . 2 ) 

where : 

= displacement of point 1 
U 2 = displacement of point 2 
Ic = stiffness coefficient. 

From [Ref. 4] the stiffness coefficient for a cylindrical 
rod with axial displacement is: 

(2.3) 

where : 

A = cross sectional area of the bar (fibers per ply) 
E = elastic modulus of fiber (axial) 

L = length of fiber. 

Combining equations (2.1a) and (2.2) yields the element 
matrix equations of: 



■ k 
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i -■ 


hi 


_-k 
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J 
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hi 



(2.4) 



for fibers, where n is the ply number. The forces and F 2 
are equal in magnitude and come from the thermal expansion of 
the material as follows. 
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The thermal strain of the material is defined as: 



e^^ = a{T-Tj (2.5) 

where : 

= thermal strain 

01 = thermal coefficient of expansion 
T = the temperature of the material 

= the temperature at which the strain is zero. 
Multiplication of equation (2.5) by Youngs modulus E yields 
the stress a, 

a =Ee=Ea {T-TJ (2.6) 

which results if free expansion or contraction is prevented. 
Multiplying equation (2.6) by the cross sectional area of the 
element yields the element forces: 

- F, =F2 =aA=AEa {T-TJ (2.7) 

at the ends of the elements where n is the ply number. 
Dropping the subscripts 1 and 2 on the forces and combining 
equation (2.7) with (2.4) and (2.3) gives the following matrix 
equations for each element: 

( 2 . 8 ) 

[-k |f|„ 

where : 

Ic = AE/L based on fiber properties and dimensions 
Ui & U 2 = displacements of points 1 & 2 for fiber n 
F = AEa(T-T„) . 
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Using the numbering system of Figure 6, the global 
stiffness matrix and global force vector are assembled as: 

(2.9) 
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where : 

kj^ = a component in the global stiffness matrix 

Ujj = a component in the global displacement vector 

Fj^ = a component in the global force vector. 

The boundary conditions imposed on the above system of 
equations are as follows. For the rigidly fixed nodal points, 

= = = = = ^ ( 2 . 10 ) 

and for the nodal points that move in unison, 

U, = U, = U, = U,= = C/2„ . (2.11) 

The system of equations of equation (2.9) with the 
boundary conditions of (2.10) and (2.11) applied were solved 
for nodal point displacements using IMSL math library 
subroutine DLSARG. 
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The stresses of the fibers per ply were calculated using 
the displacements obtained from the above system of equations 
as follows (see figure 4) : 






where 

= strain from actual motion of nodal points= (U 2 -U 3 ^) /L 
6jj, = strain causing fiber stress 

= fiber strain for free expans ion=o; (T-T„) . 
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Substitution of appropriate expressions for and gives 

an expression for mechanical strain of: 



Finally multiplying equation (2.13) by Young's modulus yields: 



where a is the stress developed in the fiber. 

Using the preceding equations a computer stress program 
was developed utilizing the Vatikiotis combustion code for 
fiber temperature profiles. The stress program determined 
fiber stresses throughout the composite material as a function 
of time and location. The stress program also determined 
failed fibers based on failure criteria of Chapter III. 
Specifics of the stress program are contained in Chapter IV. 




(2.13) 




(2.14) 
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III. FAILURE CRITERIA 



A. INTRODUCTION 

The stress analysis code determines fiber stresses 
throughout the composite plate and then determines if failure 
of any of the fibers has occurred. In this chapter possible 
fiber failure criteria are considered. The four failure modes 
considered were tensile strength failure, compressive strength 
failure, fiber buckling, and fiber consumption. Of the four 
modes considered only three were incorporated in the stress 
analysis code. 

B. TENSILE FAILURE 

The first criteria considered was failure due to exceeding 
the fibers tensile strength of 300 ksi. When subjected to 
heat graphite fibers tend to contract, which is reflected in 
their thermal expansion coefficient of - 0 . 55X10'^/°F . 
Initially, the fiber closest to the heat source will increase 
in temperature then over time the remaining fibers will heat. 
The plate fibers will continue to rise in temperature as long 
as the heat source is applied. Since thermal strain is 
proportional to temperature, the magnitude of fiber 
contraction increases with fiber temperature. Due to the 
fibers closest to the fire being hotter than those further 
from the fire, fibers closest to the fire would, if free to do 
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so, contract more than those farther away. This variation in 
tendency to contract causes induced stresses since all fibers 
are restrained to contract uniformly as seen in Figure 4. 



In order to estimate the temperature required to reach 
tensile strength, a limiting case of Equation (2.14) was used. 
Since the fiber closest to the heat source is the first to 
heat, its temperature begins to rise before any significant 
change in the other fiber temperatures. As a worst case, this 
fiber could be heated to a very high temperature while all the 
remaining fibers remained very close to their initial 
temperature. If the number of fibers were relatively high, 
this situation could be modeled as one heated fiber with fixed 
ends (i.e. not allowed to move axially). Using this model 
equation (2.14) becomes: 

S^=-Ea(T-TJ (3.1) 

or, by rearranging, the minimum temperature required for this 
fiber to reach the tensile strength is: 



Sr 

T = — + T 

-Ea 



(3.2) 



where : 

Sj. = tensile strength of the fiber (300 )csi) 

Tmin = the minimum temperature for tensile failure 
E = elastic modulus (31.0X10^ psi) 

Tjo = reference temperature (assumed to be 60°F) 
a = thermal expansion coefficient ( - 0 . 55X10‘®/°F) . 
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Solving the above equation for the minimum fiber 
temperature gives a value of approximately 17,600°F. This 
temperature is unachievable since graphites melting point is 
approximately 6000°F and, moreover, the combustion code is 
unable track temperatures over approximately 2600°F. 

Due to the above calculation, the probability of any 
fibers failing in the tensile mode is highly unlikely. 
However, a quick check for tensile failure was incorporated 
into the computer program as verification of the above 
calculations. Tensile failure was not observed during any of 
the case studies. 

C. COMPRESSIVE STRENGTH FAILURE 

The next criteria considered was purely axial compression 
of the fibers. For the graphite fibers under consideration 
their compressive strength is 260 ksi. However, this high 
value is quite deceiving. These fibers have a very small 
cross section and relatively long length, causing them to be 
susceptible to buckling failure when in compression. Stresses 
at buckling are much lower than the compressive fracture 
strength. For this reason compressive strength of the fibers 
was not used as a failure criteria. 

D. FIBER BUCKLING 

The third failure mode considered was fiber buckling. The 
graphite fibers have a very high slenderness ratio due to 
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their small cross sectional diameter and relatively long 
length. These fibers could be susceptible to Euler column 
buckling. Evaluation of this mode of failure requires a 
"typical" length for the unsupported fibers. This fiber 
length was assumed to be 10 inches, corresponding to a 
reasonable sized hole burned through the material. These 
fibers were further assumed to have fixed end conditions. 

The development of Euler Buckling formula is available in 
several texts and will not be redeveloped here. From [Ref. 5] 
the formula for Euler column buckling assuming fixed end 
conditions with constant circular cross sectional area is; 



a 



CZ 



4L^ 



(3.3) 



where : 

= critical stress for column buckling 
E = elastic modulus of the fibers 
D = the Diameter of the fibers 
L = the length of the unsupported fibers. 

Using the above formula, the calculated value of is 

approximately 19 psi. This value is very small and quite 
conservative since the model assumes no lateral support. 

As the epoxy matrix burns, some of the combustion products 
will stay behind forming a new matrix material. This new 
matrix structure consisting of residual combustion products, 
although weaker than the original epoxy matrix, renders 
lateral support to the fibers. This additional support 
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significantly increases the critical buckling strength of the 
fibers . 

Rosen [Ref. 6] , in his work on critical buckling stress of 
composite materials, develops two analytical models for 
buckling of composites, namely; the extensional mode and the 
shear mode. Of these two modes, the mode most likely for 
fiber volume fractions of approximately 20% is the shear mode. 
The shear mode formula is: 



'fcr ' 






(3.4) 



where : 

= fiber critical buckling stress 
Vf = fiber volume fraction 
Gj„ = shear modulus of the matrix material. 

A search of literature provided little information on the 
properties of epoxy combustion products for use as a matrix 
material. However, after some consideration, this problem was 
overcome by assuming that the effect of the combustion 
products, as a matrix material, could be modeled as a small 
fraction (F) of the original graphite -epoxy composite buckling 
critical stress given by equation (3.4) . This fraction (F) is 
associated with the amount of lateral support the graphite 
fibers receive from the combustion products. Using this 
assumption the following equation was used for fiber buckling 
criteria. 
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F 



(3.5) 



0<=r 



V^{1- V,) 



where is the reduced buckling critical stress and F is 
the combustion product fraction, that is, the fraction of the 
original stress due to combustion product lateral support. 
Failure of the graphite fibers due to exceeding the reduced 
buckling stress as calculated above is used by the stress 
analysis program. Different values of the combustion product 
fraction (F) are assumed to assess the effect of this fraction 
on failure. 

E. FIBER CONSUMPTION 

Fiber consumption was considered as the last failure mode. 
As the plate burns, under some conditions, the supply of 
oxygen in the interior of the plate becomes depleted. At that 
time the burning process moves to the X=0 surface where the 
air is being drawn into the plate. From then on the plate 
continues to burn at the surface of the plate slowly 
decreasing in thickness as the temperature profile across the 
plate tends to a constant value. This mode of burning is 
surface recession. 

The relatively constant temperature profile produced 
across the plate in the surface recession mode, could reduce 
the stresses seen by the fibers. If these stresses are small 
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in magnitude, the possibility of fiber failure by consumption 
is possible. Surface recession fiber consumption was 
incorporated into the stress analysis code to eliminate burnt 
away fibers from the analysis. 
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IV. COMPUTER PROGRAMS 



A. MODIFICATIONS TO THE COMBUSTION CODE 

Only superficial modifications were made to Vatikiotis' 
combustion code for this study. First, the computer code was 
brought up to Fortran 77 language by changing a few read and 
write format statements. Second, a single statement reading 
numerous input variables from the input file was broken up 
into individual read statements in order to facilitate an easy 
to read comment type input file format. This input file 
containing a short description of each input variable is 
contained in Appendix D. Third, a few output statements were 
added to facilitate building output data files for use in the 
stress analysis. The fourth and final change to this program 
was to increase the number of nodal points so that the present 
program could handle up to 101 nodal points vice the older 31 
nodal points. Increases in computer speed along with 
significant memory resources over the past few years 
facilitated this change. The amount of accuracy was 
significantly improved by using a larger number of nodal 
points, as seen in both Figure 8 and Table 1. Both display 
the same typical run of the combustion code varying only the 
number of nodal points used. Temperatures at the locations of 
X=0" (the fire surface), X=0.25" (half the medium thickness). 
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and X=0.5" (the surface opposite the fire) are compared. As 
seen, previous models of the program were not grid independent 
for the X=0" location. 



NODAL POINT EFFECT ON TEMPERATURE 

o 




NUMBER OF NODAL POINTS USED 



Figure 7: A comparison of the combustion code number of 
nodal points and calculated temperature. 
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TABLE 1: INCREMENTAL TEMPERATURE CHANGE VERSUS NODAL POINTS 



USED. 



TIME OF RUN: 60 seconds 
HEAT FLUX APPLIED: 10,000 BTU/ft^ hr 
DIFFERENTIAL PRESSURE ACROSS THE MEDIUM: 0.8 psi 
MEDIUM THICKNESS: 0.5 inches 


NUMBER 

OF 

NODAL 

POINTS 

USED 


GRAPHITE TEMPERATURE °F AT LOCATION 


X=0 inch 


X=0.25 Inch 


X=0.5 inch 


30 


822.46 


607.35 


486.48 


31 


860.57 


625.83 


486.23 


INCREMHNTAL 

CHANGE 


4.4% 


3.0% 


0.05% 


70 


1066.23 


662.77 


493.31 


71 


1068 . 11 


662 . 17 . 


493 . 14 


INCREMENTAL 

CHANGE 


0.2% 


0.02% 


0 . 03% 


100 


1099.51 


672.21 


494.42 


101 


1100.33 


672.43 


494.54 


INCREMENTAL 

CHANGE 


0 . 07% 


0.03% 


0.02% 
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B. STRESS ANALYSIS PROGRAM 



The stress program, as shown in Figure 8, obtains 
temperature information from an output data file produced by 
Vatikiotis' combustion code. Location data from the 
combustion code was presented in normalized (X/L) form since 
the equations Vatikiotis used were non-dimensionalized. In 
the stress code these locations were converted into absolute 
location points (X) . Until burn off of the fibers starts 
(i.e. prior to surface recession) the non dimensionalized 
locations (X/L) remain fixed. However, when surface recession 
starts, combustion moves to the X=0 surface and the thickness 
of the material decreases with time. As the thickness of the 
material changes the stress code continually verifies stress 
code fiber locations removing burned away fibers from the 
stress analysis. 

Vatikiotis in his combustion code formulation used a 
method of "smearing" to account for both the fiber and air 
filled void as shown in Figure 9. Denoting the porosity of 
the medium as p, we have dV^=pdV and dVg=(l-p)dV where dV is 
the total differential volume of air and graphite fiber, dV^ 
is the differential volume of air and dVg is the differential 
volume of graphite. The effect of this smearing process is to 
replace individual fibers and voids, which are extremely large 
in number, by a homogenous material. This process is utilized 
in the basic differential equations presented in Appendix A. 
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STRESS ANALYSIS CODE FLOW CHART 
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c 
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OUTPUT FILE 




IMSL 
SUBROUTINE 
DLSARG 




SUBROUTINE 
READER 



SUBROUTINE 

INTERP 



V 



SUROUTINE MATRIX 
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SUBROUTINE 
MOD2FL 
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SUBROUTINE 
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READ IN COMBUSTION 
CODE PROFILES 



0 



INTERPOLATE FOR 
FIBER TEMP 
AND DIAMETER 



I 



^orX 



REMOVE PREVIOUSLY 
FAILED FIBERS 

T~" ‘ 

BUILD GLOBAL 
MATRIX AND VACTOR 



APPLY FIXED END 
BOUNDARY 
CONDITIONS 
(TYPE I) 

" t " 

APPLY COMMON END 
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(TYPE II) 

T 

SOLVE MODIFIED 
SYSTEM OF 
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SUBROUTINE BROKEN 



0 



STRESS OUTPUT 
FILE 



( )• 




Figure 8: Flow chart of Che Stress Analysis Code 
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•SMEARING* PROCESS USED 
IN THE COMBUSTION CODE DEVELOPMENT 



PORUS 

MEDIUM 



AIR 



GRAPHITE 







-I- 



dV 



p dV 



(1-p) dV 



Wh*r*; p i( the porosity (i.e. traction of air) 



Figure 9: The process of "smearing" used in the combustion 
code . 

The stress code accommodated the idea of a homogeneous 
material and utilized temperature output by the following 
procedure (see Figure 10) . First the graphite temperatures 
along with nodal point locations are determined from the 
combustion code. These combustion code locations are compared 
to the center location of each ply. When the ply center 
location is bracketed by the combustion code locations the 
stress code performs quadratic interpolation to determine the 
graphite fiber temperature for that ply. By this method, 
accurate temperature profiles can be obtained from the 
combustion code using 101 nodal points while at the same time 
the number fibers used in the stress code may be varied. 
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Figure 10: Interpolation method from combustion code 



temperatures to stress code temperatures. 
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Elemental matrix equations are determined using 
Equation (2.8). These elemental matrices are assembled into 
a global set of equations of the general form: 
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where : 

K = the global stiffness matrix 
U = the global displacement vector 
F = the global force vector 
m = the total number of nodal points (2n) 
j and k = are arbitrary integers where l<j and k<m. 
The boundary conditions imposed on this matrix equation (4.1) , 
fit two different forms; Type I and Type II. Their 
application is performed as follows. The Type I boundary 
conditions consist of those that fit the form: 

Uj = b (4.2) 

where : 

Uj = displacement at node j , 
b = a known constant. 
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This type of boundary condition is directly comparable to the 
fixed end of the fibers (i.e. the locations of all the odd 
numbered nodal points) where b=0. Applying a single boundary 
condition, Uj=b to the equation set (4.1) above gives the 
following modified matrix equation: 



1 - j ••• k - m 
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The stress code was written to apply each fixed end boundary 
condition using the same method as (4.3) above. 

The Type II boundary conditions fit the form of: 

Uj ^ aU^ = c (4.4) 

where : 

Uj and U)^ = are two degrees of freedom 
a = known constant 
c = known constant. 
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Application of the general boundary condition of Equation (4 .4) 
to the matrix Equation (4.1) above yields: 



1 
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The detail for developing this modified matrix equation is 
given by Akin [Ref. 8] . 

The Type II boundary condition of equation (4.2) above is 
directly relatable to the fiber ends that move in unison, 

U 2 + (-1) = 0 

+ (-1) t/g = 0 

C/g + (-1) C/g = 0 (4.6) 

U2n-2^ (- 1 ) = 0 

with a=-l and c=0. The Type II boundary conditions of 
Equation (4.6) were applied in the stress analysis program. 

The final system of equations with both sets of boundary 
conditions applied was solved for the displacement vector U 
using the IMSL library subroutine DLSARG. 

Once nodal point displacements Uj ( j =1, 2 . . . , 2n) were Icnown 
fiber stresses a were determined using equation (2.14). This 
process was continued for each time step of output data 
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obtained from the combustion code. The fiber stress analysis 
code is included in Appendix B. 
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V. CASE STUDIES 



The combuscion code requires input data to start the 
analysis. Several input data files were produced encompassing 
various combinations of composite thickness, heat flux, 
differential pressure and combustion product fraction (F) . 
Specific values of these variable parameters as used in the 
case studies are listed in Table 2. Several other important 
parameters were treated as constants for these studies as 
listed in Table 3. The input data files required by the 
combustion code were built covering all 1008 combinations of 
the variable parameters. These data files were processed by 
the combustion code and then the temperature profile data was 
further processed for stress failure information. Most of the 
lower heat flux results, primarily below 4000 btu/ hr-ft^, 
produced relatively minor increases in plate temperatures 
which resulted in insignificant stresses and no failure. But 
the majority of runs gave interesting results. The final 
results of this effort are provided in graphical form in the 
following case studies. 
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TABLE 2 : THE PARAMETERS VARIED IN THE CASE STUDIES . 



INPUT PARAMETERS THAT WERE VARIED FOR THE CASE STUDIES 
WITH THE VALUES INVESTIGATED 



COMPOSITE THICKNESS; 0.25", 0.5", 0.75", 1.0" 
CONSTANT HEAT FLUX 

(thousand btu/hr- f t^) : 1 , 1.25, 1.5, 1.75, 2, 2.5, 



3, 3. 

^ C 


5, 


4.0, 


4.5, 


5, 

C- 




^ / 


i . ^ f 


o , 8 


• / 


DIFFERENTIAL PRESSURE (psi) 


: 


0.05, 


0.5, 


0.7 



COMBUSTION PRODUCT FRACTION (%): 1, 0.75, 0.5, 0.25 



TABLE 3: PARAMETERS TREATED AS CONSTANT FOR THESE STUDIES. 

CONSTANT INPUT PARAMETERS TO TEE COMBUSTION CODE 
FOR THE CASE STUDIES 

AMBIENT AIR TEMPERATURE: 80°F 
AMBIENT AIR PRESSURE: 2116.8 Ib/ft^ 

TORTUOSITY: 1.4 

PARTICLE SHAPE FACTOR: .91 

ORDER OF THE CHEMICAL REACTION: 0.5 

TIME STEP 

(between combustion code profiles): 0.5 seconds 
STOICHIOMETRIC RATIO: 0.375 
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A. TEMPERATURE PROFILES 



As Stated in previously, induced stresses are dependent 
upon the fiber temperatures. Figures 11 through 16 show the 
temperature profiles of some representative combustion code 
runs up to the point of initial fiber failure. Note that the 
maximum temperature is always at the bottom of the plate where 
the heat flux is being applied and the air is coming in, while 
the cooler fibers are at the opposite surface. The main 
concept to bare in mind for these studies is that graphite 
fibers contract upon heating, so the hotter fibers will be in 
a tensile mode (see Figure 4) . Since the tensile strength of 
these hot fibers is very high, these fibers will not fail. 
However, the fibers that are cooler will be in a compression 
mode with a low critical buckling strength and they will fail 
by buckling. 

A comparison of Figures 11, 12, and 13 show the effect of 
increasing the heat flux on a constant thickness material. 
Notice that for the lower heat fluxes, the temperatures are 
lower and the profile is more gradual. This, of course, is 
due to the lower heat flux applied but also due to the ability 
of the air flow to convectively transfer the heat away from 
the heated surface. If the heat flux is low, the flow of air 
through the medium has more time to distribute the heat 
through the material before failure occurs. If the heat flux 
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is high, the surface temperature rises quickly and heating is 
confined primarily to the area close to the heated surface. 

A comparison of Figures 11 and 16 indicates the changes 
in the temperature profile when the highest and lowest 
differential pressures are applied for a given heat flux and 
material thickness. As stated previously, the convective heat 
transfer in the medium has a large effect on the temperature 
profiles for the lower heat fluxes. For a constant thickness, 
the pressure gradient AP/L will be higher for a higher 
differential pressure. This higher pressure gradient causes 
a higher air flow through the material which increases the 
convective heat transfer. This increased heat transfer causes 
both a lower and more gradual temperature profile in the low 
heat flux cases. However, in the high heat flux cases this 
effect is insignificant as seen in Figures 12 and 14 for a 
heat flux of 10,000 btu/hr-ft^. 

As a final comparison Figures 12 and 15 are compared. 
These figures show the dependence of the temperature profile 
on thickness of the material while holding differential 
pressure constant . The temperatures at the surface of the 
1.0" plate are higher than the 0.5" plate. This can be 
explained as follows. If the AP is a fixed value, then the 
thicker plate has a smaller pressure gradient (AP/L) because 
the thickness L is larger. This smaller AP/L gradient will 
cause less air flow through the thicker material resulting in 
a lower amount of convective heat transfer. The end result for 
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the thicker material is a higher and more steep temperature 



profile . 
compressive 



This profile causes more fibers to be in 
mode and subsequently fail by buckling. 



the 
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TEMPERATURE PROFILE 

WITH RESPECT TO TIME 



o 



TIME TO FIRST FAILURE=12.1 sec 
HEAT FLUX=10,000 BTU/HR FT^ 
0.5 INCH THICKNESS 
AP=0.7PSI 




Figure 12: Temperature profile up to first failure for 0.5 
10,000 btu/hr ft^ and 0.7 psi. 
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TEMPERATURE PRORLE 

WITH RESPECT TO TIME 






TIME TO FIRST FAILURE=12.0 sec 
HEAT FLUX=10,000 BTU/HR FT^ 
0.5 INCH THICKNESS 
AP=0.05PSI 




Figure 14: Temperature profile up to first failure for 0.5", 
10,000 btu/hr ft^ and 0.05 psi. 
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B. TIME TO INITIAL FAILURE (AP EFFECT) 

Figures 17 through 20 present the time to initial failure 
showing its dependency on both heat flux and differential 
pressure. As heat flux decreases, the time to initial failure 
progressively increases. For each specific AP there appears 
to be a lower bound of heat flux below which failure will 
never occur. The figures also show that this lower bound 
value increases with increasing AP. However, this lower bound 
cannot be exactly determined through use of a computer code. 
Through the use of finite increments of heat flux, given 
previously in Table 2, this lower bound for failure was 
roughly determined. As seen in Figures 17 through 20 this 
lower bound of heat flux was found to within 500 btu/hr-ft^. 

For all thicknesses, the effect of pressure is relatively 
small except at the lower heat fluxes where the effect of air 
flow through the material has a large effect on the convective 
heat transfer. These pressure effects become ineffectual at 
higher heat fluxes. 
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• The time uu xuxuxctj. xctxxuxi: 
differential pressure) for 1.0", and F=0.5%. 



of 




I I I I I I I I I I I I I I I I I I I — I 

0 1 2345 6789 10 

HEAT FLUX 10^ BTU/HR FT^ 



Figure 18: The time to initial failure (effect of 

differential pressure) for 0.75", and F=0.5%. 
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Ficfure 19: The time to initial failure 

differential pressure) for 0.5", and F=0.5%. 
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Figure 20: The time to initial failure (effect of 

differential pressure) for 0.25", and F=0.5%. 
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C. TIME TO INITIAL FAILURE (F EFFECT) 

Figures 21 through 24 demonstrate how the critical 
buckling strength, and hence time to initial failure, varies 
with the combustion product factor (F) . The value of F 
greatly effects not only the time to failure but determines 
whether or not failure will occur. As F increases there is 
more lateral support for the fibers so the critical buckling 
strength increases. This increase results in fewer fibers 
failing by buckling. In the case of the 0.25" plate none of 
the fibers reached the critical buckling strength when F=1.0%. 
All figures show that for a given heat flux, as F increases, 
the time to initial failure increases. 

The result of this analysis shows that the value assumed 
for the combustion product fraction (F) greatly influences the 
amount of time to initial failure of the material. The 
figures show that when F increases the resistance to buckling 
increases and buckling is less likely to occur. 
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Figure 21: The time to initial failure (effect of combustion 
product fraction) for 1.0", and AP=0.7 psi. 




Figure 22: The time to initial failure (effect of combustion 
product fraction) for 0.75", and AP=0.7 psi. 
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Figure 23: The time to initial failure (effect of combustion 
product fraction) for 0.5", and AP=0.7 psi. 




Figrure 24: The time to initial failure (effect of combustion 
product fraction) for 0.25", and AP=0.7 psi. 



D. PROGRESSION OF MATERIAL FAILURE (AP EFFECT) 

Figures 25 and 26 track the percentage of failed fibers 
from the time of the first observed fiber failures. Due to 
the 0.5 second time step used for temperature profile output, 
some of these curves start at time zero with a significant 
number failed fibers. Actually the failure represented by 
these curves commenced some time between the time of zero of 
the figures and the previous 0.5 second interval. Since 
observation of trends in this case are important, the 
increased computer time to determine a more precise starting 
time was not justified. 

The general trend is, the higher the differential 
pressure, the lower the percentage of failure. This is 
apparently is due to the more gradual temperature profiles 
caused by the higher air flows through the material. For the 
higher differential pressures/higher air flowrates, more 
fibers are heated to approximately the same temperature, and 
are therefore in the tensile mode prior to the time of initial 
fiber. These tensile fibers do not fail. 

Another item worth mentioning is the rapidity of failure. 
The materials progressed to their maximum percentage of fiber 
failure very rapidly. For the 1.0" thicker material this 
maximum failure percentage is reached within 3 seconds of the 
first fiber failure and for the 0.5" thinner material 
approximately 4 seconds. These times are independent of 
differential pressures. 
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TIME FROM FIRST FAILURE (sec) 



Figure 25: The progression of material failure (the effect 
of AP) for 0.5", and 5,000 btu/ hr-ft^. 




TIME FROM FIRST FAILURE (sec) 



Figure 26: The progression of material failure (the effect 
of AP) for 1.0", and 5,000 btu/hr-ft^. 
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E. PROGRESSION OF MATERIAL FAILURE (EFFECT OF HEAT FLUX) 

Figures 27 and 28 indicate that as heat flux increases, 
time to maximum failure changes only marginally, effecting the 
3 to 4 second failure time by only 1 second. 
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TIME FROM FIRST FAILURE (sec) 



Figure 27: The progression of material failure (the effect 
of heat flux) for 0.5", and AP=0.7 psi. 




Figure 28: The progression of material failure (the effect 
of heat flux) for 1.0", and AP=0.7. 
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F. MAXIMUM FAILURE (EFFECT OF AP) 

Figures 29 through 32 are a consolidation the maximum 
fiber failure for material failure. These maximum percentages 
varied slightly with pressure at lower heat fluxes. As shown, 
for a given material thickness, higher heat fluxes give higher 
percentage of failed fibers. This is due to the more non- 
uniform temperature profile discussed previously. The maximum 
percentage failure seen in these studies is the case of the 
one inch plate at 25,000 btu/hr-ft^, which reached a maximum 
of 97.5%. 
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Figure 29: The maximum percentage of failure (the effect of 
AP) for 1.0", and F=0.5%. 




Figure 30: The maximum percentage of failure (the effect of 
AP) for 0.75", and F=0.5%. 
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Figure 31: The maximum percentage of failure (the effect of 
AP) for 0.5", and F=0.5%. 




Figure 32: The maximum percentage of failure (the effect of 
AP) for 0.25", and F=0.5%. 
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G. MAXIMUM FAILURE (EFFECT OF F) 

The effect of the combustion product fraction is shown in 
Figures 33 and 34. As expected, assuming a higher value of 
the combustion product fraction (F) reduces the amount of 
failure because of the increase in buckling strength. For the 
0.25" plate there was no failure above a factor of 1.0% for 
range of heat fluxes investigated. 

These trends were independent of differential pressure. 
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Figure 33: The maximum percentage of failure (the effect of 
combustion product fraction) 0.5", and AP=0.7 psi. 




Figure 34: The maximum percentage of failure (the effect of 
combustion product fraction) 0.25" and AP=0.7. 
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H. MAXIMUM FAILURE (EFFECT OF THICKNESS) 

As seen in Figures 35 and 36, the thickness of the 
material has minimal effect on the final percentage of 
material failure for higher heat fluxes. However at lower 
heat fluxes the thinner materials experience significantly 
less failure. In all cases the final percentage of fiber 
failure was between 42% and 97.5%. 
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Figure 35: The maximum percentage of failure (the effect of 
material thickness) F=0.5%, AP=0.7 psi. 




Figure 36: The maximum percentage of failure (the effect of 
material thickness) F=0.5%, and AP=0.05. 
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VI. CONCLUSIONS 



The conclusions drawn from this study are as follows: 

• The amount of time required before onset of fiber failure 
ranged from a very small 3 seconds (for a 0.25" plate at 
25,000 btu/hr-ft^) up to a maximum of 210 seconds (for a 
1.0" plate at 2,000 btu/hr-ft^) . Given the reaction times 
of highly trained fire fighting teams onboard aircraft 
carriers, it is quite conceivable that failure may be 
averted for the longer onset of failure times. 

• Once the failure of carbon graphite -epoxy composite 
material begins maximum failure occurs very quickly in 
less than 4 seconds. This time is of such short duration 
that reducing the amount of maximum failure once failure 
begins is not possible. 

• The predominant factors effecting survivability of carbon- 
graphite epoxy composite materials are the heat flux 
applied, the combustion product factor and the thickness 
of the material. The differential pressure effect varies 
from very important at low heat fluxes to insignificant at 
high heat fluxes. However, the effect of differential 
pressure is mitigated by the fact that response times are 
in all cases very small at the high heat fluxes. Once 
failure occurs, maximum failures in the range of 80-90% 
preclude any future use of the material. 
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APPENDIX A; COMBUSTION CODE EQUATIONS 



LIST OF EQUATION SYMBOLS 

fj. - dynamic viscosity of air 
4> - oxygen concentration 

- air density 

Pg - graphite density 

- specific heat capacity of air 

Cg - specific heat capacity of graphite 
Dg - Mass diffusion coefficient 
h - convection heat transfer coefficient 
kg - air thermal conductivity 

kg - effective thermal conductivity of porus graphite 

kj, - pseudo thermal conductivity to account for radiation 

P - pressure 

p - porosity 

Rq 2 - oxygen reaction rate 

Rg - graphite reaction rate 

Tg - air temperature (°R) 

Tg - graphite temperature (°R) 
t - time 

u - pore velocity (of air) 

z - specific internal area 
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FINAL FIELD EQUATIONS 

Following are the equations used in the Vatikiotis 
Combustion code. These equations are presented here in their 
form prior to non-dimensionalization. These equations were 
non-dimensionalized prior to use in the actual code. 

THE GRAPHITE HEAT TRANSFER EQUATION 

J^[(l-P) iK * - 1-.) * V = 



THE INTERNAL AIR HEAT TRANSFER EQUATION 





- PPaSa^ 



dx 



hz{T^-T^) 



+ u 



d(pP) 

dx 



P a ^a 



dt 



OXYGEN TRANSPORT EQUATION 







THE INTERNAL FLOW EQUATION 

d^P ^ ^ _ H ^(PPa) 

dx^ Pa dx w dx \i dx dx p^w dt 
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BOUNDARY CONDITIONS IMPOSED 

For the X=0 location, 

<a «» 



4 > = 4 >„ 



P = 



During heat addition (1-p) {k^ + k^.) = -g. 



*st 



dT„ 



After fire is removed (k^ + k.) = hAT„ - Tj + oeiT 

® ^ dx ^ ^ 



For the X=L location. 



(k, * k,)^ ^ - TJ 



- ae{T* - T^) 



dT. 



a _ 



dx 



= 0 



dx 



= 0 



P = P„ + A P 



where : 

qg^. = the heat flux applied 

= ambient temperature (°R) 

= ambient oxygen concentration 
hj^ and h 2 = convection heat transfer coefficients 
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INITIAL CONDITIONS 



There are a number of ways to start these "combustion" 
problems. One way is to set the oxygen concentrations and the 
air and graphite temperatures to ambient conditions then apply 
an external heat source. This "natural" starting scheme 
requires some amount of computer effort to bring the plate up 
to conditions where temperatures have increased significantly. 
Instead of starting the problem at ambient conditions, 
computational effort can be reduced by starting the problem 
with the plate at elevated temperatures. 

In this particular study the initial conditions were; 

Graphite temperature 400°F 

Air temperature 400°F 

Oxygen concentration 0.0172 Vo/tt? (ambient level). 
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APPENDIX B: STRESS ANALYSIS CODE 



♦EDWARD FAXLANGER JULY 1992 

* THIS PROGRAM CALCULATES STRESSES INDUCED IN LONGITUDINAL FIBERS 

* RESULTING FROM TEMPERATURE VARIATION ACROSS THE MEDIUM 
******************************************************************** 

C DESCRIPTION OF VARIABLES 
C U = DISPLACEMENTS OF NODAL POINTS 

C BC = VALUE OF THE TYPE I BOUNDARY CONDITION FOR NODAL POINT IBC 

C A = MATRIX CONTAINING INDIVIDUAL STIFFNESS 

C B = MATRIX CONTAINING INDIVIDUAL FORCE VECTOR 

C AA = MATRIX CONTAINING GLOBAL STIFFNESS 

C BB = MATRIX CONTAINING GLOBAL FORCE VECTOR 

C Y = LOCATIONS OF THE INDIVIDUAL FIBER (FROM BOTTOM OF PLATE) 

C TEMP = INTERPOLATED TEMPERATURE OF THE FIBER 
C DIAM = INTERPOLATED DIAMETER OF THE FIBER 
C SIGMA = CALCULATED STRESS OF THE INDIVIDUAL FIBER 
C YIN = LOCATIONS OF NODAL POINTS FROM THE COMBUSTION CODE 
C TEMPIN=TEMPERATURES FROM THE COMBUSTION CODE 

C DIAMIN=THE DIAMETER AT THE NODAL POINT FROM THE COMBUSTION CODE 
C IBC = THE NUMBER OF THE NODE HAVING A TYPE I BOUNDARY CONDITION 
C NODES= A MATRIX USED FOR ELEMENT CONNECTIVITY TO THEIR NODAL POINTS 
C Y1 = LOCATION OF THE NODAL POINTS BROUGHT IN FROM VATIKIOTIS CODE 
C NFLAG = A FLAG TO INDICATE TYPE OF FAILURE 

C NELEM =A MATRIX CONTAINING THE NAMES OF THE UNBROKEN FIBERS 
C NE = NUMBER OF FIBER ELEMENTS 
C NNODE = TOTAL NUMBER OF FIBER NODES 

C PERCNT= PERCENTAGE OF ORIGINAL BUCKLING STRENGTH USED 
C TSTREN= TENSILE STRENGTH OF FIBERS 

C CSTREN= COMPRESSIVE STRENGTH OF FIBERS (NO BUCKLING) 

C EFIB = ELASTIC MODULUS OF THE FIBERS 
C VF = FIBER VOLUME FRACTION (ASSUMED CONSTANT) 

C EMAT = ELASTIC MODULUS OF THE MATRIX 
C POISM = POISON RATIO OF THE MATRIX 
C ALPHA = THERMAL EXPANSION COEFFICIENT (AXIAL) 

C WIDTH = THE ASSUMED LENGTH OF EACH FIBER (ASSUMED TO BE ONE) 

C TEMREF= REFERENCE TEMPERATURE (WHERE STRAIN EQUAL ZERO) 

C NNP = THE NUMBER OF NODAL POINTS READ IN FROM COMBUSTION CODE 

C TIME = TIME OF THE TEMPERATURE PROFILE BEING ANALYZED 
C THICK = THE ORIGINAL THICKNESS OF THE MATERIAL (BEFORE BURNING) 

C HEIGHT= THE THICKNESS OF THE MATERIAL AT THE PRESENT TIME 
NEUSED= THE NUMBER OF UNBROKEN FIBERS AT PRESENT TIME 
NELEM = A STORAGE OF UNBROKEN FIBERS NAMES 
GMAT = SHEAR MODULUS OF THE EPOXY MATRIX MATERIAL 

MAIN PROGRAM 

IMPLICIT REAL*8 (A-H,0-Z) 

C THE FOLLOWING PARAMETER N1 SETS THE DIMENSION OF ARRAYS AND 
C NEEDS TO BE EQUAL TO THE HIGHER OF 1)THE NUMBER OF NODES BEING 
C READ IN **OR** 2) THE NUMBER OF ELEMENTS BEING ANALYZED 
PARAMETER (Nl=10l) 

PARAMETER (NN=2*N1) 

C COMMON STATEMENT FOR MATERIAL PROPERTIES 
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COMMON/ PROPER/TSTREN, CSTREN, EFIB , EMAT, ALPHA, WIDTH , THICK, POISM, VF 
C 

REAL*8 U(NN) ,BC(NN) ,A(2,2),B(2) ,BB(NN) ,AA(NN,NN) ,Y(N1) , 
&TEMP(N1) ,DIAM(N1) , SIGMA (Nl) ,YIN(N1) ,TEMPIN(N1) ,DIAMIN(N1) 
DIMENSION IBC(NN) , NODES (NN, NN) ,INDEX(2) ,NFLAG(N1) ,NELEM(N1) 

*THE FOLLOWING 3 STATEMENTS INCREASE THE NUMBER OF EQUATIONS 
*THAT THE EQUATION SOLVER DLSARG WILL HANDLE. IF LARGER SPACE IS 

* NEEDED DLSARG WILL GIVE AN ERROR MESSAGE AND RECOMMEND LARGER 

* WORK SPACES 

COMMON/ WORKS P /RWKS P 
REAL RWKSP (82236) 

CALL IWKIN(82236) 

C DATA FILE FOR OUTPUT OF DATA 

OPEN(UNIT=ll, FILE=' /TEST31 STRESS B1 ' , STATUS= ' OLD ' ) 

C DATA FILE FOR READING INPUT DATA 

OPEN (UNIT=12 , FILE=' /TEST31 TEMPER B1 ' , STATUS= ' OLD ' ) 

C NE IS THE NUMBER OF FIBERS ASSUMED FOR ANALYSIS 
NE = 30 

NNODE=NE*2 

C ENTER THE PERCENTAGE OF MATRIX SUPPORT 
PERCNT=.005 

C FIBER CHARACTERISTICS; LONGITUDINAL COMPRESSIVE AND TENSILE 
C STRENGTHS (PSD , ELASTIC MODULUS AND THERMAL COEFFICIENT 
TSTREN=300 .0D3 
CSTREN=-260.0D3 
EFIB=31.0D6 
EMAT= . 5D6 
VF= .196 
POISM=.35 
ALPHA=-0.55D-6 

C WIDTH OF SIMULATED HOLE (LENGTH OF FIBERS) ASSUMED TO BE UNITY 
C SINCE IT HAS NO EFFECT ON PROBLEM 
WIDTH=1 . ODO 

C REFERENCE TEMPERATURE IN DEGREES F WHERE INITIAL STRESSES WERE ZERO 
TEMREF= 60. DO 

C INITIALIZE FLAGS TO ZERO (NFLAG KEEPS TRACK OF BROKEN FIBERS) 

C (NSTOP STOPS THE PROGRAM ) 

NSTOP=0 
DO M=1,NE 
NFLAG (M) =0 
END DO 



* READ IN OF INFORMATION FROM COMBUSTION CODE OUTPUT 
READ(12,*)NNP 
CALL CHECK (NNP, NSTOP) 

I F (NSTOP. NE.O) GOTO 100 

CALL READER (Nl , NNP, TIME , HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
IF (NSTOP. NE.O) GOTO 100 

*LOOP FOR FILLING ALL ELEMENT HEIGHTS ASSUMING EVEN SPACING 
THICK=HEIGHT 
YINC=THICK/ (NE) 

YSTART=YINC/2 
DO M=1,NE 

Y(M) =YSTART+YINC* (M-1) 

END DO 
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200 BURN=THICK-HEIGHT 
C INITIALIZE ARRAYS TO ZERO 
DO 1=1, NN 
BB (I) =0 . ODO 
DO J=1,NN 
AA(I, J) =0 .ODO 
END DO 
END DO 

DO 1=1, NE 

C IF THE FIBERS ARE ALREADY CONSUMED BY COMBUSTION (OUTSIDE EDGE) 

C NULL THEIR EFFECT BY SETTING THEIR FLAG TO 1 
IF(Y(I) ,LT.BURN)THEN 
NFLAG(I) =1 
END IF 

C DETERMINE THE TEMPERATURE AND DIAMETER OF NON- FAILED FIBERS THIS 
C LOOP IS FOR INTERPOLATION OF DATA RECEIVED FROM "STRESS INPUT" FILE. 
C THE RECEIVED DATA IS IN RELATIVE SPACING AND MUST BE CHANGED TO 
C ABSOLUTE SPACING IN ORDER TO DETERMINE PERCENTAGE FAILURE OF FIBERS 
C INTERPOLATE NON-BROKEN FIBERS ONLY 
IF(NFLAGd) .EQ.O) THEN 
DO N=l,NNP-2 

IF( (Y(I) .GE.YIN(N) +BURN) .AND. (Y(I) . LE . YIN (N+2 ) +BURN) ) THEN 
CALL INTERP(YIN(N) +BURN, YIN (N+1 ) +BURN, YIN (N+2 ) +BURN, 

&DIAMIN(N) ,DIAMIN(N+1) ,DIAMIN(N+2) ,Y(I) ,DIAM(I) ) 

CALL INTERP (YIN (N) +BURN, YIN (N+1 ) +BURN, YIN (N+2 ) +BURN, 

&TEMPIN(N) ,TEMPIN(N+1) ,TEMPIN(N+2) ,Y(I) ,TEMP(I) ) 

END IF 
END DO 
END IF 
END DO 

C LOOP FOR DETERMINING THE ELEMENTS THAT ARE NOT BROKEN (IE THE 
C ONES THAT ARE TO BE ANALIZED) 

NEUSED=0 
DO 1=1, NE 

IF (NFLAG (I) .EQ. 0) THEN 
NEUSED=NEUSED+1 
NELEM(NEUSED) =I 
END IF 
END DO 

IF (NEUSED . LT . 2 ) THEN 

WRITE (11, *) 'ALL BUT TWO ELEMENTS HAVE BROKEN' 

WRITE (11, *) 'ANALYSIS IS STOPPED/PROGRAM STOPPED' 

GOTO 100 
END IF 

NNUSED=NEUSED*2 

* LOOP FOR ESTABLISHING CONNECTIVITY OF NODES TO THEIR ELEMENTS 

* FOR USE IN BUILDING GLOBAL MATRIX AND GLOBAL VECTOR 

L=1 

DO 1=1, NEUSED 
DO J=l,2 
NODES (I, J) =L 
L=L+1 
END DO 
END DO 



* LOOP TO BUILD GLOBAL MATRIX (AA) AND VECTOR (BB) 
DO IE=1, NEUSED 
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CALL MATRIX (A, B, TEMP (NELEM( IE) ) -TEMREF, DIAM (NELEM ( IE) ) ) 

INDEX (1) =NODES (IE, 1) 

INDEX ( 2 ) =NODES ( IE , 2 ) 

DO 1=1,2 
II=INDEX(I) 

BB (ID =BB (II) +B (I) 

DO J=l,2 
JJ=INDEX( J) 

AA(II, JJ) =AA(II, JJ) +A(I, J) 

END DO 
END DO 
END DO 

♦LOOP FOR TYPE I BOUNDARY CONDITION INITIALIZATION 
DO I=1,NNUSED-1, 2 
IBC (I) =I 
BC(I) =0.0D0 
END DO 

* APPLY TYPE 1 BOUNDARY CONDITIONS 
DO I=1,NNUSED-1,2 
J=IBC (I) 

DO K=1,NNUSED 

BB(K) =BB (K) -BC (I) *AA(K, J) 

AA(K, J)=0.0 
AA(J,K) =0.0 
AA(J, J) =1.0 
BB (J) =BC (I) 

END DO 
END DO 

C APPLY TYPE 2 BOUNDARY CONDITIONS 
DO J=2,NNUSED-2, 2 

CALL MOD2FL(NN,NNUSED, J, J+2, - 1 . ODO , 0 . ODO , AA, BB) 

END DO 

C SOLVE THE SYSTEM OF EQUATIONS USING IMSL LIBRARY 
CALL DLSARG(NNUSED,AA, NN,BB,1,U) 

DO 1=1, NE 
SIGMA(I) =0 .ODO 
END DO 

DO I=1,NEUSED 
STRAIN=U (1*2) /WIDTH 

SIGMA (NELEM (I) ) =EFIB* (STRAIN-ALPHA* (TEMP (NELEM(I) ) -TEMREF) ) 

END DO 

C LOOP TO CALCULATE STRESSES OF EACH ELEMENT AND OUTPUT RESULTS 
C TO OUTPUT FILE 

WRITE (11, 9) TIME 
WRITE (11,6) 

DO 1=1, NE 

IF(NFLAGd) .NE.0)WRITE(11,8) I,NFLAG(I) 

IF(NFLAGd) .EQ.0)WRITE(11,7) I,SIGMA(I) ,DIAM(I) 

END DO 

6 FORMAT (5X, 'FIBER NUMBER' , 5X, 'FIBER STRESS (PSI )', 2X, ' DIAMETER' ) 

7 FORMAT(12X,I3,11X,G10.4,6X,G10.4,3X,G11.4,2X,I1) 

8 FORMAT (17X, 13, 4X, 'PREVIOUSLY BROKEN, NFLAG=',I1) 

9 FORMATdX, 'TIME=' ,F8.3, ' SECONDS') 

CALL BROKEN(Nl,NE, SIGMA, DIAM, NFLAG, BROKE, PERCNT) 

WRITE (11, *) '*************i*f***************************' 
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WRITE (11, *) 

CALL READER (N1 , NNP, TIME , HEIGHT, YIN, TEMP IN, DIAMIN, NSTOP) 
IF(NSTOP.NE.O)GOTO 100 
GOTO 200 
100 END 



SUBROUTINE MATRIX (A, B , TEMP , DIAM) 

THIS SUBROUTINE CALCULATES THE ELEMENT LEVEL 
STIFFNESS MATRIX (A) AND VECTOR (B) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/PROPER/TSTREN,CSTREN, EFIB , EMAT, ALPHA, WIDTH, THICK, POISM, VF 
REAL*8 A(2,2) ,B(2) 

PI=4*ATAN(1 .0) 

AREA=DIAM* DIAM* P 1/4 . 0 

S=AREA*EFIB/WIDTH 

A(l, 1) =S 

A(l,2) =-S 

A(2 , 1) =-S 

A(2,2) =S 

B (2) =EFIB*AREA*ALPHA*TEMP 
B(l) =-B(2) 

END 



SUBROUTINE BROKEN (N1 , NE , STRESS , DIAM, NFLAG , BROKE , PERCNT) 

THIS SUBROUTINE DETERMINES THE NUMBER OF BROKEN FIBERS (BROKE) 

DUE TO TENSILE AND COMPOSITE BUCKLING 
NFLAG EQUALS 0 : FOR UNFAILED FIBER 

l:FOR BURNT AWAY FIBERS 
2: FOR OVER STRESS BY TENSION 
3 : FOR OVER STRESS BY BUCKLING 
IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/ PROPER/TSTREN , CSTREN, EFIB , EMAT, ALPHA, WIDTH, THICK, POISM, VF 
REAL*8 STRESS (Nl) , DIAM (Nl) 

DIMENSION NFLAG (Nl) 

PI=4*ATAN(1.0) 

C CALCULATE COMPOSITE BUCKLING STRENGTH OF THE FIBER 
C AS A FRACTION OF THE ORIGINAL COMPOSITE BUCKLING STRENGTH 
GMAT=EMAT/ (2 . * (1 . +POISM) ) 

BUCKLE=-GMAT/ (1 . -VF) *PERCNT 
DO 1=1, NE 

C DO NOT CHECK PREVIOUSLY BROKEN FIBERS 
IF(NFLAGd) .EQ.O) THEN 
C CHECK FOR TENSILE STRESS FAILURE 

IF (STRESS (I) .GE.TSTREN) THEN 
NFLAG (I) =2 
END IF 

C CHECK FOR BUCKLING FAILURE 

IF (STRESS (I) .LE. BUCKLE) THEN 
NFLAG (I) =3 
END IF 
END IF 
END DO 

C SUMMATION OF EACH TYPE OF BREAKAGE 
BURNT=0 . 0 
TENSIL=0.0 
BUCK=0 . 0 
DO 1=1, NE 

IF(NFLAGd) .EQ.l) BURNT=BURNT+1 
IF(NFLAGd) .EQ.2) TENSIL=TENSIL+1 
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IF(NFLAGd) .EQ.3) BUCK=BUCK+1 
END DO 

TBROKE =BURNT+TENS I L+BUCK 
WRITE (11,5) BURNT*100/NE 
5 FORMAT (IX, ' PERCENTAGE BURNT FIBERS' , F5 . 1) 

WRITE (11,4)TENSIL*100/NE 

4 FORMAT (IX, ' PERCENTAGE BROKEN IN TENSILE MODE',F5.1) 
WRITE (11 , 2) BUCK*100/NE 

2 FORMAT (IX, 'PERCENTAGE OF BUCKLED FIBERS ', F5 . 1 ) 

WRITE (11, 1) TBROKE*! 00/NE 

1 FORMAT (IX, 'TOTAL PERCENTAGE OF FAILED FIBERS ', F5 . 1 ) 
END 



SUBROUTINE INTERP (XI , X2 , X3 , Y1 , Y2 , Y3 , XX, YY) 

THIS SUBROUTINE TAKES 3 X VALUES, 3 Y VALUES. THEN USING A GIVE XX 
VALUE PERFORMS QUADRATIC INTERPOLATION TO DETERMINE THE YY VALUE. 
IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 A(3,3) ,B(3) ,C(3) 

A(1,1)=1.D0 
A(2, 1) =1 .DO 
A(3, 1) =1 .DO 
A(l, 2) =X1 
A(2,2) =X2 
A(3,2) =X3 
A(l, 3) =X1**2 . 

A(2,3) =X2**2 . 

A(3, 3) =X3**2 . 

B (1) =Y1 
B (2) =Y2 
B (3) =Y3 

CALLING IMSL LIBRARY SUBROUTINE TO SOLVE A (3 , 3 ) *C ( 3 ) =B (3 ) 

TO DETERMINE THE QUADRATIC COEFFICIENTS 
CALL DLSARG(3,A,3,B,1,C) 

YY=C (1) +C(2) *XX+C(3) *XX**2 . 

END 

SUBROUTINE MOD2FL (NN, NSIZE , LI , L2 , A, B , S , C) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION S(NN,NN), C (NN) 

APPLY TYPE 2 CONSTRAINTS TO FULL EQUATIONS 
METHOD AND SUBROUTINE BY: J. ED AKIN 
************************** 

NN = THE TOTAL SIZE OF THE PASSED 
NSIZE = WORKING SIZE OF ARRAYS 
S = FULL SQUARE MATRIX 
C = COLUMN MATRIX 

L1,L2 = CONSTRAINED DEGREES OF FREEDOM SUBCRIPTS OF D MATIRX 
A,B = CONSTRAINT COEFFICIENTS 

D(L1) + A*D(L2) = B , S*D = C 

INITIAL CALCULATIONS 
S12 = S(L1,L2) 

Cll = C(L1) 

SUBTRACT A* COLUMN LI FROM COLUMN L2 

SUBTRACT B* COLUMN LI FROM C 

DO 10 I = 1, NSIZE 

S(I,L2) = S(I,L2) - A*S(I,L1) 

S(L2,I) = S(I,L2) 

C(I) = C(I) - B*S(I,L1) 
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S (I,L1) = 0.0 
S(L1,I) = 0.0 

10 CONTINUE 

S(L1,L1) =1.0 

S(L2,L2) = S(L2,L2) - A*S12 + A*A 

S(L2,L1) = A 

S(L1,L2) = A 
C(L1) = B 

C(L2) = C(L2) - A*C11 + A*B 

RETURN 
END 
C 
C 

SUBROUTINE READER (N1 , NNP , TIME , HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
C THIS SUBROUTINE READS IN THE NODAL POINT VALUES FROM THE 
C COMBUSTION CODE 

IMPLICIT REAL*8 (A-H,0-2) 

REAL*8 YIN(Nl) , TEMPIN (Nl) , DIAMIN (Nl) 

READ(12,*,END=10) 

READ(12,*,END=10) 

READ(12,*,END=10) 

READ (12, *,END=10) TIME 
READ(12,*,END=10) 

READ(12,*,END=10) HEIGHT 
READ (12 , * , END=10) 

C CONVERT HEIGHT TO INCHES 
HEIGHT=HEIGHT*12 . 0 
DO 1=1, NNP 

READ(12, *,END=10) K,Y1N(I) ,TEMPIN(I) ,DIAMIN(I) 

C CONVERT INPUT LOCATION (YIN) AND DIAMETER (DIAMIN) TO INCHES 
YIN(I) =YIN(I) *12 . 0 
DIAMIN(I) =DIAMIN(I) *12.0 
END DO 
GOTO 20 
10 NSTOP =1 
2 0 END 
C 
C 

SUBROUTINE CHECK (NNP , NSTOP) 

C QUICK CHECK TO SEE IF THE NUMBER OF INPUT NODAL POINTS IS ODD 
C SO QUADRATIC INTERPOLATION CAN BE USED. IF NOT STOP PROGRAM 
MCHECK=MOD (NNP , 2 ) 

IF (MCHECK.EQ. 0) THEN 

WRITE (*, 2) 'THE NUMBER OF INPUT NODAL POINTS MUST BE AN ODD' 
WRITE(*,2) 'NUMBER IN ORDER FOR THIS PROGRAM TO USE QUADRATIC' 
WRITE (*, 2) 'INTERPOLATION. CHECK NNP IN INPUT DATA' 

WRITE (*, 2) 'THIS PROGRAM IS STOPPED' 

2 FORMAT (IX, A) 

NSTOP =1 
END IF 
END 
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APPENDIX C: SAMPLE OUTPUT DATA PROM STRESS CODE 



TIME= 0.000 SECONDS 
FIBER NUMBER 



FIBER STRESS (PSI) 



1 


0 .OOOOE+00 


2 


0 .OOOOE + 00 


3 


0 .OOOOE + 00 


4 


O.OOOOE+00 


5 


O.OOOOE+00 


6 


0 .OOOOE+00 


7 


O.OOOOE+00 


8 


0 .OOOOE+00 


9 


0 .OOOOE + 00 


10 


0 . OOOOE+00 


11 


O.OOOOE+00 


12 


O.OOOOE+00 


13 


0 . OOOOE+00 


14 


0 .OOOOE+00 


15 


0 .OOOOE+00 


16 


0 .OOOOE+00 


17 


0 .OOOOE + 00 


18 


0 .OOOOE+00 


19 


0 .OOOOE+00 


20 


0 .OOOOE+00 


21 


0 .OOOOE+00 


22 


0 .OOOOE+00 


23 


0 .OOOOE+00 


24 


0 . OOOOE+00 


25 


0 .OOOOE+00 


26 


O.OOOOE+00 


27 


0 .OOOOE+00 


28 


0 .OOOOE+00 


29 


0 .OOOOE+00 


30 


0 .OOOOE+00 


BURNT FIBERS 


0.0 



0.0 



PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 0.0 
TOTAL PERCENTAGE OF FAILED FIBERS 0 . 0 



DIAMETER 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0 .8333E-02 
0 . 8333E-02 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 . 8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0 .8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0 .8333E-02 
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TIME= 8.810 SECONDS 
FIBER NUMBER 



FIBER STRESS (PSI) 



1 


2777 . 


2 


2236. 


3 


1791. 


4 


1385 . 


5 


1022 . 


6 


705.4 


7 


433.8 


8 


206.0 


9 


18.97 


10 


-131.2 


11 


-249.3 


12 


-340.1 


13 


-408.6 


14 


-459.2 


15 


-495.9 


16 


-522.0 


17 


-540.2 


18 


-552.6 


19 


-561.0 


20 


-566.6 


21 


-570.3 


22 


-572 . 6 


23 


-574 . 1 


24 


-575 . 1 


25 


-575.7 


26 


-576.0 


27 


-576.3 


28 


-576.4 


29 


-576.5 


30 


-576.5 


BURNT FIBERS 


0.0 



PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 0.0 
TOTAL PERCENTAGE OF FAILED FIBERS 



0.0 



0.0 



★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 



DIAMETER 
0 .8333E-02 
0 . 8333E-02 
0 .8333E-02 
0 . 8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0 . 8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0 .8333E-02 
0.8333E-02 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
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TIME= 17.810 SECONDS 
FIBER NUMBER 
1 
2 

3 

4 

5 

6 

7 

8 
9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

PERCENTAGE BURNT FIBERS 
PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 13.3 
TOTAL PERCENTAGE OF FAILED FIBERS 13.3 



FIBER STRESS (PSI) 
3524 . 

3022 . 

2597. 

2190. 

1805. 

1443 . 

1107 . 

797 . 8 

515.9 
261.7 
34.97 

-165.1 

-339.8 

-490.6 

-619.4 

-728.2 

-819.1 

-894.3 

-955.8 

-1006. 

-1046 . 

-1077 . 

- 1102 . 

-1122 . 

-1136 . 

-1148 . 

-1156 . 

-1161 . 

-1165. 

-1167 . 

0.0 

0.0 






DIAMETER 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0 . 8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0 .8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0 . 8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 

0.8333E-02 
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TIME= 26.560 SECONDS 



FIBER NUMBER 




FIBER STRESS 


(PSD 


1 




1389 . 




2 




1519 . 




3 




1520. 




4 




1392 . 




5 




1174 . 




6 




902.0 




7 




603.2 




8 




297.3 




9 




-3.231 




10 




-290.8 




11 




-561.4 




12 




-812.8 




13 




-1044 . 




14 




-1255. 




15 




-1445 . 




16 




-1616. 




17 




-1768 . 






18 


PREVIOUSLY 


BROKEN 




19 


PREVIOUSLY 


BROKEN 




20 


PREVIOUSLY 


BROKEN 




21 


PREVIOUSLY 


BROKEN 




22 


PREVIOUSLY 


BROKEN 




23 


PREVIOUSLY 


BROKEN 




24 


PREVIOUSLY 


BROKEN 




25 


PREVIOUSLY 


BROKEN 




26 


PREVIOUSLY 


BROKEN 




27 


PREVIOUSLY 


BROKEN 




28 


PREVIOUSLY 


BROKEN 




29 


PREVIOUSLY 


BROKEN 




30 


PREVIOUSLY 


BROKEN 



PERCENTAGE BURNT FIBERS 0.0 
PERCENTAGE BROKEN IN TENSILE MODE 0.0 
PERCENTAGE OF BUCKLED FIBERS 56.7 
TOTAL PERCENTAGE OF FAILED FIBERS 56.7 



DIAMETER 
0.8333E-02 
0 .8333E-02 
0 .8333E-02 
0 .8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0.8333E-02 
0 .8333E-02 
0.8333E-02 
0 .8333E-02 
0.8333E-02 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
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APPENDIX D; SAMPLE INPUT DATA TO COMBUSTION CODE 



80. EO TAMB AMBIENT TEMPERATURE 

50.4E0 UINF DIFF. PRESS. ACROSS MEDIUM LB/FT2 OR SURF VELOCITY 

2116. 8E0 PAMB AMBIENT PRESSURE LB/FT2 

1.4E0 TORT TORTUOSITY OF THE MEDIUM 



0 . 868055555555555507E-03 DINT 
0 .868055555555555507E-03 SINIT 
0 .416666666666666667E-01 XLENGTH 
0 .416666666666666667E-01 XINIT 



LUMP FIBER DIAMETER PER CELL 
DISTANCE BETWEEN LUMPED FIBERS 
TOTAL MEDIUM THICKNESS 
INITIAL THICKNESS OF THE MEDIUM 



53 .34E0 


RGAS 


86.0E0 


GCOND 


0 .231E0 


SPHTG 


70.3E0 


ROG 


0.9E0 


EMIS 


.91E0 


PHI 


l.OEO 


HLENTH 


14085. EO 


ENTHAL 


0 .5E0 


ORDNTH 


0.5E0 


EXTEMP 


0 .375E0 


STIOCH 


2.065E6 


RECOEFF 


28836 .EO 


ENERGY 


30000 .EO 


SQ 


4 .DO 


T 


O.EO 


TSTART 


200 .EO 


TSTOP 


40 .EO 


TQ 


l.OOD-5 


HSTART 


1 .E-10 


HMIN 


0.500 


HMAX 


0.0010 


RMSEPS 


2500 .EO 


TMAX 


50. EO 


TMIN 


.020E0 


CMAX 


- .OOIEO 


CMIN 


O.EO 


SLl 


O.EO 


SL2 


O.EO 


SL3 



6 

0 

1 

25 

0 

0 

0 

0 

0 

12 

1 

0 

500 

0 



MAXDER 

IPRT 

ITIME 

NNP 

lELEM 

IDIAM 

I SPACE 

lUNIF 

I RE AC 

IBC 

lie 

IQ 

NSTOP 

ITRIP 



THE GAS CONSTANT FOR AIR 
THERMAL CONDUCTIVITY OF GRAPHITE 
SPECIFIC HEAT OF GRAPHITE 
DENSITY OF GRAPHITE 
THE EMISSIVITY OF THE GRAPHITE 
SHAPE FACTOR 1.0 FOR SPHERE .91 FOR CYLINDERS 
REF LENGTH FOR SURFACE HEAT TRANSFER (NORMALLY 1) 
THE ENTHALPY OF CHEMICAL REACTION 
THE ORDER OF THE CHEMICAL REACTION 
THE EXPONENT OF THE TEMPERATURE FOR REACTION RATE 
STOICHIOMETRIC RATIO FOR THE REACTION 
ARRHENIUS EQUATION COEFFICIENT 

ACTIVATION ENERGY DIVIDED BY THE IDEAL GAS CONSTANT 
THE HEAT FLUX (SIMULATING THE FIRE) 

THE APPROXIMATE TIME INTERVAL FOR PRINTING RESULTS 
THE STARTING TIME OF THE PROBLEM 
THE STOP TIME OF THE PROBLEM 
THE HEATER TURN OFF TIME (FIRE IS OUT) 

THE INITIAL TIME INTERVAL IN SUBROUTINE SDESOL 
THE MINIMUM TIME INTERVAL FOR SUBROUTINE SDESOL 
THE MAXIMUM TIME INTERVAL FOR SUBROUTINE SDESOL 
ERROR CRITERIA FOR SUBROUTINE SDESOL 
THE MAX TEMP ON PRINTER PLOTS (Y AXIS) 

THE MIN TEMP ON PRINTER PLOTS (Y AXIS) 

THE MAX OXYGEN CONC . ON PRINTER PLOTS (Y AXIS) 

THE MIN OXYGEN CONC. ON PRINTER PLOTS (YAXIS) 

SLOPE FOR A LINEAR PROFILE FOR GRAPHITE TEMP 
SLOPE FOR A LINEAR PROFILE FOR AIR TEMP 
SLOPE FOR A LINEAR PROFILE FOR CONCENTRATION 
MAX ORDER EQUATION SOLVER USED IN SDESOL 
PRINT PARAMETER FOR SDESOL 

MULTIPLIER OF PRINT TIME STEP ONCE PROGRAM IS IN SURF RECESS 
NUMBER OF NODAL POINTS 

0=AUTO 1=USER INPUT (SEE BELOW) OF NODAL LOCATIONS 
0=AUTO 1=USER INPUT (SEE BELOW) OF PARTICLE DIAMETER 

0=UINF IS DIFFERENTIAL PRESSURE 1=UNIF IS SURFACE VELOCITY 
0=AUTO STOICH. RATIO OF CARBON OXIDES 1=USER ENTRY OF RATIO 
BOUNDARY CONDITION PARAMETER 

l=AUTO LINEAR TEMPG, TEMPA, CONC. 02 0=USER ENTRY (SEE BELOW) 
LOCATION OF HEATER 0=HEATER AT ZERO 1=HEATER AT XLENGTH 
MAX NUMBER OF TIMES ALLOWED THROUGH INTEGRATION LOOPS 
LEAVE ZERO UNLESS RESTARTING PROGRAM IN SURF RECESSION MODE 



800. EO --GRAPHITE TEMP NODE 1] *THESE 3 ITEMS ARE REQUIRED ONLY 
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800 .EO 
. 0172E0 
0 .DO 
.OOIDO 
.005D0 
. OIDO 
. 025D0 
.05D0 
. IDO 
. 15D0 
.2D0 
.25D0 
.3D0 
.35D0 
.4D0 
.45D0 
.5D0 
.55D0 
. 6D0 
.65D0 
. 7D0 
. 738D0 
. 8D0 
. 85D0 
.9D0 
.95D0 
l.DO ** 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800. DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800. DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800. DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 
800 .DO 



- AIR TEMP NODE 1 ] *IF IIC=1 . THEY MUST FOLLOW "ITRIP" OR 

--OXYGEN CONC NODE 1 - - ] *MANUAL ENTRY OF NODE LOCATIONS IF USED 
******THE FOLLOWING DATA IS THE X LOCATION OF THE NODES 
*THESE LOCATIONS ARE ONLY NEEDED IF IELEM=1 
^REQUIRING NODAL LOCATION INPUT. IF IELEM=0 
*THE COMPUTER CALCULATES A LINEAR SPACING OF THE 
*NODES . 

* 

* 

* IF THIS MANUAL ENTRY IS USED IT MUST 

* FOLLOW INPUT OF THE VARIABLE "ITRIP" 

* AND HAVE A VALUE FOR EACH NODAL POINT 



*****Ein3 OF NODAL POINT LOCATION INPUT*************** 

$$$$$$$$ START OF MANUAL ENTRY OF GRAPHITE TEMP FOR EACH NODE 
$ 

$$$$$THIS INFORMATION ON TEMPERATURE AND CONCENTRATION 
IS ONLY REQUIRED IF IIC=0 REQUIRING MANUAL DATA INPUT 
FOR EACH NODAL POINT. THE $$$$$$ DATA MUST FOLLOW DATA 
FOR "ITRIP" IF NODAL LOCATIONS ARE NOT BEING MANUALLY 
INPUT (IF IELEM=l)OR ELSE THIS DATA FOLLOWS NODAL POINT 
MANUAL ENTRY LOCATIONS 



DATA MUST BE GIVEN FOR EACH NODAL POINT 



$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$$$$$$LAST MANUAL ENTRY OF GRAPHITE TEMPERATURE%%%%%%%% 
$$$$$$START MANUAL ENTRY OF AIR TEMPERATURE EACH NODE## 
$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 
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800. DO 

800. DO 

800 .DO 

800 .DO 

800 .DO 

800. DO 

800 .DO 

800. DO 

800 .DO 

800 .DO 

800 .DO 

800. DO 

800 .DO 

800 .DO 

800 .DO 

800. DO 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

-0075D0 

. 0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

. 0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

.0075D0 

. 0075D0 



$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$$$$$LAST MANUAL ENTRY OF AIR TEMPERATURE EACH NODE 
$$$$START MANUAL ENTRY OF OXYGEN CONCENTRATION EACH NODE 
$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$$$$$LAST MANUAL ENTRY OF OXYGEN CONCENTRATION EACH NODE 
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APPENDIX E: SAMPLE OUTPUT DATA FROM COMBUSTION CODE 



31 NODAL POINTS 



AT TIME= 

0.00 S 

MEDIUM THICKNESS= 

4166666666667D-01 . 4166666666667D-01 
NP X TEMP-G 



1 O.OOOOOOOOOOOOOOOOE+00 

2 0 .1388888888888889E-02 

3 0 .2777777777777778E-02 

4 0 .4166666666666666E-02 

5 0 -5555555555555555E-02 

6 0 . 6944444444444443E-02 

7 0.8333333333333331E-02 

8 0 .9722222222222221E-02 

9 0 . llllllllllllllllE-01 

10 0 . 1250000000000000E-01 

11 0 . 1388888888888889E-01 

12 0 . 1527777777777777E-01 

13 0.1666666666666666E-01 

14 0 .1805555555555555E-01 

15 0 . 1944444444444444E-01 

16 0 .2083333333333333E-01 

17 0 .2222222222222222E-01 

18 0 . 2361111111111111E-01 

19 0.2499999999999999E-01 

20 0.2638888888888888E-01 

21 0 . 2777777777777777E-01 

22 0.2916666666666666E-01 

23 0.3055555555555555E-01 

24 0 .3194444444444444E-01 

25 0.3333333333333333E-01 

26 0.3472222222222222E-01 

27 0 .3611111111111110E-01 

28 0.3749999999999999E-01 

29 0.3888888888888888E-01 

30 0.4027777777777777E-01 

31 0.4166666666666666E-01 



0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000 000 000 0000 OOE+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0.4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0.4000000000000000E+03 
0 .4000 000 000 000 OOOE+03 
0 .400 0000 000 000 OOOE+03 
0 .4000000000000000E+03 
0 .400 0000 000 0000 OOE+03 
0 .400 0000 000 000 OOOE+03 
0 .4000000000000000E+03 
0 .4000000000000000E+03 
0 .4000 000 000 000 OOOE+03 
0 .4000 000 000 0000 OOE+03 
0 .4000 000 000 000 OOOE+03 
0 .400 0000 00 0000 OOOE+03 
0 . 4000000000000000E+03 
0 .4000 000 000 0000 OOE+03 
0.4000000000000000E+03 
0 .4000 0000 00 0000 OOE+03 



DIAMETER 

0 . 69444444444 44444E- 03 
0 .6944444444444444E-03 
0 . 6944444444444444E-03 
0 .6944444444444444E-03 
0 .6944444444444444E-03 
0 . 6944444444444444E-03 
0 .6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944 44 444444444 4E- 03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 .694444 444 4444 444 E- 03 
0 . 6944444444444444E-03 
0 .69444 4444 444 4444E- 03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 694444 4444444444E- 03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
0 . 6944444444444444E-03 
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AT TIME= 

8.81 S 

MEDIUM THICKNESS = 

4166666666667D-01 .41666 6666 6667D- 01 
NP X TEMP-G 



1 0 .OOOOOOOOOOOOOOOOE+00 

2 0.1388888888888889E-02 

3 0 . 2777777777777778E-02 

4 0 .4166666666666666E-02 

5 0.5555555555555555E-02 

6 0 .6944444444444443E-02 

7 0.8333333333333331E-02 

8 0 .9722222222222221E-02 

9 O.llllllllllllllllE-01 

10 0.1250000000000000E-01 

11 0.1388888888888889E-01 

12 0 .1527777777777777E-01 

13 0 .1666666666666666E-01 

14 0 .1805555555555555E-01 

15 0 .1944444444444444E-01 

16 0 .2083333333333333E-01 

17 0 .2222222222222222E-01 

18 0 .2361111111111111E-01 

19 0 .2499999999999999E-01 

20 0 .2638888888888888E-01 

21 0 .2777777777777777E-01 

22 0 .2916666666666666E-01 

23 0 .3055555555555555E-01 

24 0 .3194444444444444E-01 

25 0 .3333333333333333E-01 

26 0 .3472222222222222E-01 

27 0 .3611111111111110E-01 

28 0 .3749999999999999E-01 

29 0.3888888888888888E-01 

30 0 .4027777777777777E-01 

31 0 .4166666666666666E-01 



0 .617165541474 3030E+03 
0 .57819 84237740606E+03 
0 .5511844 753523755E+03 
0. 52 613082 74 8824 58E+03 
0 .5035640136879362E+03 
0.4836232366654954E+03 
0 .4663593 795839235E+03 
0 . 451724909249414 8E+03 
0. 43958334 609 1013 8E+03 
0.4297262006976060E+03 
0.4218950356130471E+03 
0.4 158052 0689 13094E+03 
0 .4111679694792949E+03 
0 .4 077084937514 894E+03 
0 .4 0517854 13371006E+03 
0 .4033636888937181E+03 
0 .4020858170994085E+03 
0 .401202019 8152163E+03 
0 .4006011841982683E+03 
0. 4001993 564763142E+03 
0 .3999347556308410E+03 
0 .3997630187488629E+03 
0 .3996530135102078E+03 
0 .3995833611737365E+03 
0 .3995396822558629E+03 
0 .3995124994579232E+03 
0 .3994956957918823E+03 
0 .3994854173856393E+03 
0 .3994793187448693E+03 
0 .3994760668482828E+03 
0 . 3994750372395511E+03 



DIAMETER 

0 . 6944444 02 154 0727E- 03 
0 .6944444306055943E-03 
0 . 6944444380346299E-03 
0 .69444444134 89135E- 03 
0. 694444442 842569 7E- 03 
0 .6944444435425437E-03 
0 . 6944444438859773E-03 
0 . 694444444 0630632E- 03 
0 . 6944444441 590080E- 03 
0 . 69444444421342 59E- 03 
0 . 6944444442455372E-03 
0 . 694444444265 104 OE- 03 
0 .6944444442773239E-03 
0 . 6944444442850919E-03 
0 . 69444444429009 14E- 03 
0 . 6944444442933 369E- 03 
0 . 69444444429 54 588E- 03 
0 .69444444429 68574E- 03 
0 . 694444444297790 6E- 03 
0 .6944444442984252E-03 
0 .69444444429 88702E- 03 
0 . 6944444442991956E-03 
0 .6944444442994469E-03 
0 .6944444442996526E-03 
0 .694444444299 8301E- 03 
0 .6944444442999913E-03 
0 .69444444430014 14E- 03 
0 .6944444443002882E-03 
0 . 6944444443004240E-03 
0 .6944444443005773E-03 
0 .6944444443006692E-03 
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AT TIME= 

17.81 S 

MEDIUM THICKNESS= 

4166666666667D-01 . 4 1666666 66667D- 01 
NP X TEMP-G 



1 0 .OOOOOOOOOOOOOOOOE+00 

2 0.1388888888888889E-02 

3 0.2777777777777778E-02 

4 0.4166666666666666E-02 

5 0.5555555555555555E-02 

6 0 . 6944444444444443E-02 

7 0.8333333333333331E-02 

8 0 .9722222222222221E-02 

9 0 . llllllllllllllllE-01 

10 0 .1250000000000000E-01 

11 0.1388888888888889E-01 

12 0 . 1527777777777777E-01 

13 0.1666666666666666E-01 

14 0.1805555555555555E-01 

15 0 . 1944444444444444E-01 

16 0.2083333333333333E-01 

17 0 . 2222222222222222E-01 

18 0.2361111111111111E-01 

19 0.2499999999999999E-01 

20 0 .2638888888888888E-01 

21 0 .2777777777777777E-01 

22 0 .2916666666666666E-01 

23 0 .3055555555555555E-01 

24 0 .3194444444444444E-01 

25 0.3333333333333333E-01 

26 0 .3472222222222222E-01 

27 0 .3611111111111110E-01 

28 0.3749999999999999E-01 

29 0.3888888888888888E-01 

30 0 .4027777777777777E-01 

31 0 .4166666666666666E-01 



0 . 69444832723 7619 8E+03 
0 . 65829524834 89 8 12E+03 
0 . 633 0284 155199 859E+03 
0 . 6085590421656422E+03 
0 .5853 068 87244 08 80E+03 
0 .56339 15122869859E+03 
0 .5429230303669187E+03 
0 .5239878470762538E+03 
0 .5066449982116862E+03 
0 .4909237200244683E+03 
0 .47682263 85980613E+03 
0 .4643106220267744E+03 
0 .4533292103391948E+03 
0 .443 79 63 8472023 81E+03 
0 .4356113 128753 018E+03 
0 .42865963 6512 99 55E+03 
0 .4228188585575783E+03 
0 .4179634359788205E+03 
0 .4139692729 89 6838E+03 
0 . 410717418276353 OE+03 
0.4 0809 68793508916E+03 
0 .4060065623 12 6001E+03 
0 .4043564178987307E+03 
0 .4030679223461711E+03 
0.4020740463757375E+03 
0 .4013 188723 800529E+03 
0 .4 00757013 7582 123E+03 
0 .4003529789248906E+03 
0. 400080599 6527675E+03 
0 .39992265943 25125E+03 
0 . 399 870699 855263 6E+03 



DIAMETER 

0 . 694444 02 0 005179 7E- 03 
0 . 6944442750792861E-03 
0. 694444354 73 35164E-03 
0. 69444439 6918603 8E-03 
0 . 6944444188805062E-03 
0 .6944444303859835E-03 
0 . 6944444 3 64 794 849E- 03 
0 . 6944444397590113E-03 
0 . 6944444415612797E-03 
0 . 69444444 25767153E- 03 
0. 694444443 165029 8E- 03 
0 .6944444435161306E-03 
0 . 69444444 3 73 2 039 8E- 03 
0 . 694444443 8 6872 54E- 03 
0 . 694444443 9 57 63 09E- 03 
0 . 6944444440168827E-03 
0 . 6944444440572177E-03 
0 . 6944444440851723E-03 
0 . 6944444 44 104 83 60E- 03 
0 . 6944444441188342E-03 
0 . 6944444441288955E-03 
0 . 69444444413 6183 IE -03 
0 . 694444444 14 14 949 E- 03 
0 . 6944444441453872E-03 
0 .6944444441482521E-03 
0 .6944444441503701E-03 
0 .694444444 1519 3 71E- 03 
0 .694444444 153 0992E- 03 
0 .694444444 153 9341E- 03 
0 .6944444441545535E-03 
0 .6944444441548335E-03 
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AT TIME= 

26.56 S 

MEDIUM THICKNESS= 

4166666666667D-01 . 4166666666667D- 01 
NP X TEMP-G 



1 O.OOOOOOOOOOOOOOOOE+00 

2 0.1388888888888889E-02 

3 0 . 2777777777777778E-02 

4 0 . 4 166666666666666E - 02 

5 0 .5555555555555555E-02 

6 0 . 6944444444444443E-02 

7 0 -8333333333333331E-02 

8 0 -9722222222222221E-02 

9 O.llllllllllllllllE-01 

10 0 .1250000000000000E-01 

11 0 .1388888888888889E-01 

12 0.1527777777777777E-01 

13 0.1666666666666666E-01 

14 0 .1805555555555555E-01 

15 0 . 1944444444444444E-01 

16 0.2083333333333333E-01 

17 0 .2222222222222222E-01 

18 0 . 2361111111111111E-01 

19 0 .2499999999999999E-01 

20 0 .2638888888888888E-01 

21 0 .2777777777777777E-01 

22 0 .2916666666666666E-01 

23 0 .3055555555555555E-01 

24 0 .3194444444444444E-01 

25 0.3333333333333333E-01 

26 0 .3472222222222222E-01 

27 0 .3611111111111110E-01 

28 0 .3749999999999999E-01 

29 0 .3888888888888888E-01 

30 0 .4027777777777777E-01 

31 0 .4166666666666666E-01 



0 .6306595038438272E+03 
0 . 639924132064 1830E+03 
0 . 6447192511167083E+03 
0 .6406985836760575E+03 
0 .6302556153620338E+03 
0 . 6156188718308147E+03 
0 .5986698701861527E+03 
0 .580797904 0862 021E+03 
0 .5629290453016817E+03 
0 .5456325613672719E+03 
0 .5292359938247957E+03 
0 . 5139161110240054E+03 
0 .4997594668844450E+03 
0 .48679850024 14993E+03 
0 .4750316490868472E+03 
0 .4644343493352758E+03 
0 .4549 653630104 84 lE+03 
0 .44 6570932 6513 69 7E+03 
0 .4391880102692778E+03 
0 .432 74710849 03999E+03 
0 .427174967 1005763E+03 
0 .42239 70689333552E+03 
0 .41833998166534 12E+03 
0 .4149 33494362 344 6E+03 
0 .4121125315457292E+03 
0 .4098188488450426E+03 
0 .4080025338059548E+03 
0 .4066233572360078E+03 
0 .40565200655722 66E+03 
0 .4050713541567405E+03 
0 .4048 77367 13 7059 7E+03 



DIAMETER 

0 . 6944429267898032E-03 
0 . 6944437360540249E-03 
0 .6944440181040645E-03 
0 .6944441916810881E-03 
0 . 694444295779 19 66E- 03 
0 . 6944443573613218E-03 
0 .69444439 3301502 IE -03 
0 .6944444 14 112 7005E- 03 
0 .6944444261547878E-03 
0 . 6944444331631678E-03 
0 . 6944444 3 72 8 7 164 5E- 03 
0 . 6944444 39 75 07 13 OE- 03 
0 . 6944444412492834E-03 
0 .69444444 2 179 59 3 IE -03 
0 . 69444444 2 7 69 82 8 7E- 03 
0 . 6944444431527831E-03 
0 . 6944444434068504E-03 
0 . 6944444435790780E-03 
0 . 6944444436982122E-03 
0 . 6944444437821565E-03 
0 . 6944444438422836E-03 
0 . 6944444438859 626E- 03 
0 . 69444444 39 18 0605E- 03 
0 .69444444 39418518E- 03 
0 .6944444439595736E-03 
0 .6944444439727761E-03 
0 .6944444439825334E-03 
0 .6944444439895999E-03 
0 .6944444439944511E-03 
0 . 6944444439974719E-03 
0 .69444444 399 857 87E- 03 
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